home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 22 / Cream of the Crop 22.iso / math / ast53src.zip / CALC.C < prev    next >
C/C++ Source or Header  |  1996-09-29  |  30KB  |  901 lines
  1. /*
  2. ** Astrolog (Version 5.30) File: calc.c
  3. **
  4. ** IMPORTANT NOTICE: The graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1996 by Walter D. Pullen
  6. ** (Astara@msn.com, http://www.magitech.com/~cruiser1/astrolog.htm).
  7. ** Permission is granted to freely use and distribute these routines
  8. ** provided one doesn't sell, restrict, or profit from them in any way.
  9. ** Modification is allowed provided these notices remain with any
  10. ** altered or edited versions of the program.
  11. **
  12. ** The main planetary calculation routines used in this program have
  13. ** been Copyrighted and the core of this program is basically a
  14. ** conversion to C of the routines created by James Neely as listed in
  15. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  16. ** available from Matrix Software. The copyright gives us permission to
  17. ** use the routines for personal use but not to sell them or profit from
  18. ** them in any way.
  19. **
  20. ** The PostScript code within the core graphics routines are programmed
  21. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  22. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  23. **
  24. ** The extended accurate ephemeris databases and formulas are from the
  25. ** calculation routines in the program "Placalc" and are programmed and
  26. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  27. ** (alois@azur.ch). The use of that source code is subject to
  28. ** regulations made by Astrodienst Zurich, and the code is not in the
  29. ** public domain. This copyright notice must not be changed or removed
  30. ** by any user of this program.
  31. **
  32. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  33. ** X Window graphics initially programmed 10/23-29/1991.
  34. ** PostScript graphics initially programmed 11/29-30/1992.
  35. ** Last code change made 9/22/1996.
  36. */
  37.  
  38. #include "astrolog.h"
  39.  
  40.  
  41. /*
  42. ******************************************************************************
  43. ** House Cusp Calculations.
  44. ******************************************************************************
  45. */
  46.  
  47.  
  48. /* This is a subprocedure of ComputeInHouses(). Given a zodiac position,  */
  49. /* return which of the twelve houses it falls in. Remember that a special */
  50. /* check has to be done for the house that spans 0 degrees Aries.         */
  51.  
  52. int HousePlaceIn(rDeg)
  53. real rDeg;
  54. {
  55.   int i = 0;
  56.  
  57.   rDeg = Mod(rDeg + 0.5/60.0/60.0);
  58.   do {
  59.     i++;
  60.   } while (!(i >= cSign ||
  61.       (rDeg >= chouse[i] && rDeg < chouse[Mod12(i+1)]) ||
  62.       (chouse[i] > chouse[Mod12(i+1)] &&
  63.       (rDeg >= chouse[i] || rDeg < chouse[Mod12(i+1)]))));
  64.   return i;
  65. }
  66.  
  67.  
  68. /* For each object in the chart, determine what house it belongs in. */
  69.  
  70. void ComputeInHouses()
  71. {
  72.   int i;
  73.  
  74.   for (i = 0; i <= cObj; i++)
  75.     inhouse[i] = HousePlaceIn(planet[i]);
  76. }
  77.  
  78.  
  79. /* This house system is just like the Equal system except that we start */
  80. /* our 12 equal segments from the Midheaven instead of the Ascendant.   */
  81.  
  82. void HouseEqualMidheaven()
  83. {
  84.   int i;
  85.  
  86.   for (i = 1; i <= cSign; i++)
  87.     chouse[i] = Mod(is.MC-270.0+30.0*(real)(i-1));
  88. }
  89.  
  90.  
  91. /* Compute the cusp positions using the Alcabitius house system. */
  92.  
  93. void HouseAlcabitius()
  94. {
  95.   real rDecl, rSda, rSna, r;
  96.   int i;
  97.  
  98.   rDecl = RAsin(RSin(is.OB) * RSinD(is.Asc));
  99.   r = -RTan(AA) * RTan(rDecl);
  100.   rSda = DFromR(RAcos(r));
  101.   rSna = rDegHalf - rSda;
  102.   chouse[sLib] = DFromR(is.RA) - rSna;
  103.   chouse[sSco] = DFromR(is.RA) - rSna*2.0/3.0;
  104.   chouse[sSag] = DFromR(is.RA) - rSna/3.0;
  105.   chouse[sCap] = DFromR(is.RA);
  106.   chouse[sAqu] = DFromR(is.RA) + rSda/3.0;
  107.   chouse[sPis] = DFromR(is.RA) + rSda*2.0/3.0;
  108.   for (i = sLib; i <= sPis; i++)
  109.     chouse[i] = Mod(chouse[i]+is.rSid);
  110.   for (i = sAri; i <= sVir; i++)
  111.     chouse[i] = Mod(chouse[i+6]+rDegHalf);
  112. }
  113.  
  114.  
  115. /* This is a new house system similar in philosophy to Porphyry houses.   */
  116. /* Instead of just trisecting the difference in each quadrant, we do a    */
  117. /* smooth sinusoidal distribution of the difference around all the cusps. */
  118.  
  119. void HousePorphyryNeo()
  120. {
  121.   real delta;
  122.   int i;
  123.  
  124.   delta = (MinDistance(is.MC, is.Asc) - rDegQuad)/4.0;
  125.   chouse[sLib] = Mod(is.Asc+rDegHalf); chouse[sCap] = is.MC;
  126.   chouse[sAqu] = Mod(chouse[sCap] + 30.0 + delta   + is.rSid);
  127.   chouse[sPis] = Mod(chouse[sAqu] + 30.0 + delta*2 + is.rSid);
  128.   chouse[sSag] = Mod(chouse[sCap] - 30.0 + delta   + is.rSid);
  129.   chouse[sSco] = Mod(chouse[sSag] - 30.0 + delta*2 + is.rSid);
  130.   for (i = sAri; i < sLib; i++)
  131.     chouse[i] = Mod(chouse[i+6]-rDegHalf);
  132. }
  133.  
  134.  
  135. /* The "Whole" house system is like the Equal system with 30 degree houses, */
  136. /* where the 1st house starts at zero degrees of the sign of the Ascendant. */
  137.  
  138. void HouseWhole()
  139. {
  140.   int i;
  141.  
  142.   for (i = 1; i <= cSign; i++)
  143.     chouse[i] = Mod((SFromZ(is.Asc)-1)*30+ZFromS(i));
  144. }
  145.  
  146.  
  147. /* In "null" houses, the cusps are always fixed to start at their cor-    */
  148. /* responding sign, i.e. the 1st house is always at 0 degrees Aries, etc. */
  149.  
  150. void HouseNull()
  151. {
  152.   int i;
  153.  
  154.   for (i = 1; i <= cSign; i++)
  155.     chouse[i] = Mod(ZFromS(i));
  156. }
  157.  
  158.  
  159. /* Calculate the house cusp positions, using the specified algorithm. */
  160.  
  161. void ComputeHouses(housesystem)
  162. int housesystem;
  163. {
  164.   char sz[cchSzDef];
  165.  
  166.   if (RAbs(AA) > RFromD(rDegQuad-rAxis) && housesystem < 2) {
  167.     sprintf(sz,
  168.       "The %s system of houses is not defined at extreme latitudes.",
  169.       szSystem[housesystem]);
  170.     PrintError(sz);
  171.     Terminate(tcFatal);
  172.   }
  173.   switch (housesystem) {
  174.   case  1: HouseKoch();           break;
  175.   case  2: HouseEqual();          break;
  176.   case  3: HouseCampanus();       break;
  177.   case  4: HouseMeridian();       break;
  178.   case  5: HouseRegiomontanus();  break;
  179.   case  6: HousePorphyry();       break;
  180.   case  7: HouseMorinus();        break;
  181.   case  8: HouseTopocentric();    break;
  182.   case  9: HouseAlcabitius();     break;
  183.   case 10: HouseEqualMidheaven(); break;
  184.   case 11: HousePorphyryNeo();    break;
  185.   case 12: HouseWhole();          break;
  186.   case 13: HouseNull();           break;
  187.   default: HousePlacidus();
  188.   }
  189. }
  190.  
  191.  
  192. /*
  193. ******************************************************************************
  194. ** Star Position Calculations.
  195. ******************************************************************************
  196. */
  197.  
  198. /* This is used by the chart calculation routine to calculate the positions */
  199. /* of the fixed stars. Since the stars don't move in the sky over time,     */
  200. /* getting their positions is mostly just reading info from an array and    */
  201. /* converting it to the correct reference frame. However, we have to add    */
  202. /* in the correct precession for the tropical zodiac, and sort the final    */
  203. /* index list based on what order the stars are supposed to be printed in.  */
  204.  
  205. void ComputeStars(SD)
  206. real SD;
  207. {
  208.   int i, j;
  209.   real x, y, z;
  210.  
  211.   /* Read in star positions. */
  212.  
  213.   for (i = 1; i <= cStar; i++) {
  214.     x = rStarData[i*6-6]; y = rStarData[i*6-5]; z = rStarData[i*6-4];
  215.     planet[oNorm+i] = RFromD(x*rDegMax/24.0+y*15.0/60.0+z*0.25/60.0);
  216.     x = rStarData[i*6-3]; y = rStarData[i*6-2]; z = rStarData[i*6-1];
  217.     if (x < 0.0) {
  218.       neg(y); neg(z);
  219.     }
  220.     planetalt[oNorm+i] = RFromD(x+y/60.0+z/60.0/60.0);
  221.     /* Convert to ecliptic zodiac coordinates. */
  222.     EquToEcl(&planet[oNorm+i], &planetalt[oNorm+i]);
  223.     planet[oNorm+i] = Mod(DFromR(planet[oNorm+i])+rEpoch2000+SD);
  224.     planetalt[oNorm+i] = DFromR(planetalt[oNorm+i]);
  225.     ret[oNorm+i] = RFromD(rDegMax/26000.0/365.25);
  226.     starname[i] = i;
  227.   }
  228.  
  229.   /* Sort the index list if -Uz, -Ul, -Un, or -Ub switch in effect. */
  230.  
  231.   if (us.nStar > 1) for (i = 2; i <= cStar; i++) {
  232.     j = i-1;
  233.  
  234.     /* Compare star names for -Un switch. */
  235.  
  236.     if (us.nStar == 'n') while (j > 0 && NCompareSz(
  237.       szObjName[oNorm+starname[j]], szObjName[oNorm+starname[j+1]]) > 0) {
  238.       SwapN(starname[j], starname[j+1]);
  239.       j--;
  240.  
  241.     /* Compare star brightnesses for -Ub switch. */
  242.  
  243.     } else if (us.nStar == 'b') while (j > 0 &&
  244.       rStarBright[starname[j]] > rStarBright[starname[j+1]]) {
  245.       SwapN(starname[j], starname[j+1]);
  246.       j--;
  247.  
  248.     /* Compare star zodiac locations for -Uz switch. */
  249.  
  250.     } else if (us.nStar == 'z') while (j > 0 &&
  251.       planet[oNorm+starname[j]] > planet[oNorm+starname[j+1]]) {
  252.       SwapN(starname[j], starname[j+1]);
  253.       j--;
  254.  
  255.     /* Compare star declinations for -Ul switch. */
  256.  
  257.     } else if (us.nStar == 'l') while (j > 0 &&
  258.       planetalt[oNorm+starname[j]] < planetalt[oNorm+starname[j+1]]) {
  259.       SwapN(starname[j], starname[j+1]);
  260.       j--;
  261.     }
  262.   }
  263. }
  264.  
  265.  
  266. /*
  267. ******************************************************************************
  268. ** Chart Calculation.
  269. ******************************************************************************
  270. */
  271.  
  272. /* Given a zodiac degree, transform it into its Decan sign, where each    */
  273. /* sign is trisected into the three signs of its element. For example,    */
  274. /* 1 Aries -> 3 Aries, 10 Leo -> 0 Sagittarius, 25 Sagittarius -> 15 Leo. */
  275.  
  276. real Decan(deg)
  277. real deg;
  278. {
  279.   int sign;
  280.   real unit;
  281.  
  282.   sign = SFromZ(deg);
  283.   unit = deg - ZFromS(sign);
  284.   sign = Mod12(sign + 4*((int)RFloor(unit/10.0)));
  285.   unit = (unit - RFloor(unit/10.0)*10.0)*3.0;
  286.   return ZFromS(sign)+unit;
  287. }
  288.  
  289.  
  290. /* Transform spherical to rectangular coordinates in x, y, z. */
  291.  
  292. void SphToRec(r, azi, alt, rx, ry, rz)
  293. real r, azi, alt, *rx, *ry, *rz;
  294. {
  295.   real rT;
  296.  
  297.   *rz = r *RSinD(alt);
  298.   rT  = r *RCosD(alt);
  299.   *rx = rT*RCosD(azi);
  300.   *ry = rT*RSinD(azi);
  301. }
  302.  
  303.  
  304. #ifdef PLACALC
  305. /* Compute the positions of the planets at a certain time using the Placalc */
  306. /* accurate formulas and ephemeris. This will supersede the Matrix routine  */
  307. /* values and is only called with the -b switch is in effect. Not all       */
  308. /* objects or modes are available using this, but some additional values    */
  309. /* such as Moon and Node velocities not available without -b are. (This is  */
  310. /* the one place in Astrolog which calls the Placalc package functions.)    */
  311.  
  312. void ComputePlacalc(t)
  313. real t;
  314. {
  315.   int i;
  316.   real r1, r2, r3, r4;
  317.  
  318.   /* We can compute the positions of Sun through Pluto, Chiron, the four */
  319.   /* asteroids, Lilith, and the (true or mean) North Node using Placalc. */
  320.   /* The other objects must be done elsewhere.                           */
  321.  
  322.   for (i = oSun; i <= oLil; i++) {
  323.     if ((ignore[i] && i > oMoo) ||
  324.       (us.fPlacalcAst && FBetween(i, oCer, oVes)))
  325.       continue;
  326.     if (FPlacalcPlanet(i, t*36525.0+2415020.0, us.objCenter != oEar,
  327.       &r1, &r2, &r3, &r4)) {
  328.  
  329.       /* Note that this can't compute charts with central planets other */
  330.       /* than the Sun or Earth or relative velocities in current state. */
  331.  
  332.       planet[i]    = Mod(r1 + is.rSid);
  333.       planetalt[i] = r2;
  334.       ret[i]       = RFromD(r3);
  335.  
  336.       /* Compute x,y,z coordinates from azimuth, altitude, and distance. */
  337.  
  338.       SphToRec(r4, planet[i], planetalt[i],
  339.         &spacex[i], &spacey[i], &spacez[i]);
  340.     }
  341.   }
  342.  
  343.   /* If heliocentric, move Earth position to object slot zero. */
  344.  
  345.   i = us.objCenter != oEar;
  346.   if (i) {
  347.     planet[oEar] = planet[oSun];
  348.     planetalt[oEar] = planetalt[oSun];
  349.     ret[oEar] = ret[oSun];
  350.     spacex[oEar] = spacex[oSun];
  351.     spacey[oEar] = spacey[oSun];
  352.     spacez[oEar] = spacez[oSun];
  353.   }
  354.   spacex[i] = spacey[i] = spacez[i] =
  355.     planet[i] = planetalt[i] = ret[i] = 0.0;
  356.   if (us.objCenter < oMoo)
  357.     return;
  358.  
  359.   /* If other planet centered, shift all positions as in Matrix formulas. */
  360.  
  361.   for (i = 0; i <= oLil; i++) if (!FIgnore(i)) {
  362.     spacex[i] -= spacex[us.objCenter];
  363.     spacey[i] -= spacey[us.objCenter];
  364.     spacez[i] -= spacez[us.objCenter];
  365.     ProcessPlanet(i, 0.0);
  366.   }
  367.   spacex[us.objCenter] = spacey[us.objCenter] = spacez[us.objCenter] =
  368.     planet[us.objCenter] = planetalt[us.objCenter] = ret[us.objCenter] = 0.0;
  369. }
  370. #endif
  371.  
  372.  
  373. /* This is probably the main routine in all of Astrolog. It generates a   */
  374. /* chart, calculating the positions of all the celestial bodies and house */
  375. /* cusps, based on the current chart information, and saves them for use  */
  376. /* by any of the display routines.                                        */
  377.  
  378. real CastChart(fDate)
  379. bool fDate;
  380. {
  381.   CI ci;
  382.   real housetemp[cSign+1], Off = 0.0, vtx, j;
  383.   int i, k;
  384.  
  385.   /* Hack: Time zone +/-24 means to have the time of day be in Local Mean */
  386.   /* Time (LMT). This is done by making the time zone value reflect the   */
  387.   /* logical offset from GMT as indicated by the chart's longitude value. */
  388.  
  389.   if (RAbs(ZZ) == 24.0)
  390.     ZZ = DegToDec(DecToDeg(OO)/15.0);
  391.   ci = ciCore;
  392.  
  393.   if (MM == -1) {
  394.  
  395.     /* Hack: If month is negative, then we know chart was read in through a  */
  396.     /* -o0 position file, so the planet positions are already in the arrays. */
  397.  
  398.     is.MC = planet[oMC]; is.Asc = planet[oAsc];
  399.   } else {
  400.     for (i = 0; i <= cObj; i++) {
  401.       planet[i] = planetalt[i] = 0.0;    /* On ecliptic unless we say so.  */
  402.       ret[i] = 1.0;                      /* Direct until we say otherwise. */
  403.     }
  404.     Off = ProcessInput(fDate);
  405.     ComputeVariables(&vtx);
  406.     if (us.fGeodetic)               /* Check for -G geodetic chart. */
  407.       is.RA = RFromD(Mod(-OO));
  408.     is.MC  = CuspMidheaven();       /* Calculate our Ascendant & Midheaven. */
  409.     is.Asc = CuspAscendant();
  410.     ComputeHouses(us.nHouseSystem); /* Go calculate house cusps. */
  411.  
  412.     /* Go calculate planet, Moon, and North Node positions. */
  413.  
  414.     ComputePlanets();
  415.     if (!ignore[oMoo] || !ignore[oNod] || !ignore[oSou] || !ignore[oFor]) {
  416.       ComputeLunar(&planet[oMoo], &planetalt[oMoo],
  417.         &planet[oNod], &planetalt[oNod]);
  418.       ret[oNod] = -1.0;
  419.     }
  420.  
  421.     /* Compute more accurate ephemeris positions for certain objects. */
  422.  
  423. #ifdef PLACALC
  424.     if (us.fPlacalc)
  425.       ComputePlacalc(is.T);
  426. #endif
  427.     if (!us.fPlacalc) {
  428.       planet[oSou] = Mod(planet[oNod]+rDegHalf);
  429.       ret[oSou] = ret[oNod] = RFromD(-0.053);
  430.       ret[oMoo] = RFromD(12.5);
  431.     }
  432.  
  433.     /* Calculate position of Part of Fortune. */
  434.  
  435.     j = planet[oMoo]-planet[oSun];
  436.     if (us.nArabicNight < 0)
  437.       neg(j);
  438.     j = RAbs(j) < rDegQuad ? j : j - RSgn(j)*rDegMax;
  439.     planet[oFor] = Mod(j+is.Asc);
  440.  
  441.     /* Fill in "planet" positions corresponding to house cusps. */
  442.  
  443.     planet[oVtx] = vtx; planet[oEP] = CuspEastPoint();
  444.     for (i = 1; i <= cSign; i++)
  445.       planet[cuspLo + i - 1] = chouse[i];
  446.     if (!us.fHouseAngle) {
  447.       planet[oAsc] = is.Asc; planet[oMC] = is.MC;
  448.       planet[oDes] = Mod(is.Asc + rDegHalf);
  449.       planet[oNad] = Mod(is.MC + rDegHalf);
  450.     }
  451.     for (i = oFor; i <= cuspHi; i++)
  452.       ret[i] = RFromD(rDegMax);
  453.   }
  454.  
  455.   /* Go calculate star positions if -U switch in effect. */
  456.  
  457.   if (us.nStar)
  458.     ComputeStars(us.fSidereal ? 0.0 : -Off);
  459.  
  460.   /* Transform ecliptic to equatorial coordinates if -sr in effect. */
  461.  
  462.   if (us.fEquator)
  463.     for (i = 0; i <= cObj; i++) if (!ignore[i]) {
  464.       planet[i]    = RFromD(Tropical(planet[i]));
  465.       planetalt[i] = RFromD(planetalt[i]);
  466.       EclToEqu(&planet[i], &planetalt[i]);
  467.       planet[i]    = DFromR(planet[i]);
  468.       planetalt[i] = DFromR(planetalt[i]);
  469.     }
  470.  
  471.   /* Now, we may have to modify the base positions we calculated above */
  472.   /* based on what type of chart we are generating.                    */
  473.  
  474.   if (us.fProgress && us.fSolarArc) {  /* Are we doing -p0 solar arc chart? */
  475.     for (i = 0; i <= cObj; i++)
  476.       planet[i] = Mod(planet[i] + (is.JDp - is.JD) / us.rProgDay);
  477.     for (i = 1; i <= cSign; i++)
  478.       chouse[i]  = Mod(chouse[i]  + (is.JDp - is.JD) / us.rProgDay);
  479.     }
  480.   if (us.nHarmonic > 1)            /* Are we doing a -x harmonic chart?     */
  481.     for (i = 0; i <= cObj; i++)
  482.       planet[i] = Mod(planet[i] * (real)us.nHarmonic);
  483.   if (us.objOnAsc) {
  484.     if (us.objOnAsc > 0)           /* Is -1 put on Ascendant in effect?     */
  485.       j = planet[us.objOnAsc]-is.Asc;
  486.     else                           /* Or -2 put object on Midheaven switch? */
  487.       j = planet[-us.objOnAsc]-is.MC;
  488.     for (i = 1; i <= cSign; i++)   /* If so, rotate the houses accordingly. */
  489.       chouse[i] = Mod(chouse[i]+j);
  490.   }
  491.  
  492.   /* Check to see if we are -F forcing any objects to be particular values. */
  493.  
  494.   for (i = 0; i <= cObj; i++)
  495.     if (force[i] != 0.0) {
  496.       planet[i] = force[i]-rDegMax;
  497.       planetalt[i] = ret[i] = 0.0;
  498.     }
  499.  
  500.   ComputeInHouses();        /* Figure out what house everything falls in. */
  501.  
  502.   /* If -f domal chart switch in effect, switch planet and house positions. */
  503.  
  504.   if (us.fFlip) {
  505.     for (i = 0; i <= cObj; i++) {
  506.       k = inhouse[i];
  507.       inhouse[i] = SFromZ(planet[i]);
  508.       planet[i] = ZFromS(k)+MinDistance(chouse[k], planet[i]) /
  509.         MinDistance(chouse[k], chouse[Mod12(k+1)])*30.0;
  510.     }
  511.     for (i = 1; i <= cSign; i++) {
  512.       k = HousePlaceIn(ZFromS(i));
  513.       housetemp[i] = ZFromS(k)+MinDistance(chouse[k], ZFromS(i)) /
  514.         MinDistance(chouse[k], chouse[Mod12(k+1)])*30.0;
  515.     }
  516.     for (i = 1; i <= cSign; i++)
  517.       chouse[i] = housetemp[i];
  518.   }
  519.  
  520.   /* If -3 decan chart switch in effect, edit planet positions accordingly. */
  521.  
  522.   if (us.fDecan) {
  523.     for (i = 0; i <= cObj; i++)
  524.       planet[i] = Decan(planet[i]);
  525.     ComputeInHouses();
  526.   }
  527.  
  528.   ciCore = ci;
  529.   return is.T;
  530. }
  531.  
  532.  
  533. /* Calculate the position of each planet with respect to the Gauquelin      */
  534. /* sectors. This is used by the sector charts. Fill out the planet position */
  535. /* array where one degree means 1/10 the way across one of the 36 sectors.  */
  536.  
  537. void CastSectors()
  538. {
  539.   int source[MAXINDAY], type[MAXINDAY], occurcount, division, div,
  540.     i, j, s1, s2, ihouse, fT;
  541.   real time[MAXINDAY], rgalt1[objMax], rgalt2[objMax],
  542.     azi1, azi2, alt1, alt2, lon, lat, mc1, mc2, d, k;
  543.  
  544.   /* If the -l0 approximate sectors flag is set, we can quickly get rough   */
  545.   /* positions by having each position be the location of the planet as     */
  546.   /* mapped into Placidus houses. The -f flip houses flag does this for us. */
  547.  
  548.   if (us.fSectorApprox) {
  549.     ihouse = us.nHouseSystem; us.nHouseSystem = 0;
  550.     not(us.fFlip);
  551.     CastChart(fTrue);
  552.     not(us.fFlip);
  553.     us.nHouseSystem = ihouse;
  554.     return;
  555.   }
  556.  
  557.   /* If not approximating sectors, then they need to be computed the formal */
  558.   /* way: based on a planet's nearest rising and setting times. The code    */
  559.   /* below is similar to ChartInDayHorizon() accessed by the -Zd switch.    */
  560.  
  561.   fT = us.fSidereal; us.fSidereal = fFalse;
  562.   lon = RFromD(Mod(Lon)); lat = RFromD(Lat);
  563.   division = us.nDivision * 4;
  564.   occurcount = 0;
  565.  
  566.   /* Start scanning from 18 hours before to 18 hours after the time of the */
  567.   /* chart in question, to find the closest rising and setting times.      */
  568.  
  569.   ciCore = ciMain; ciCore.tim -= 18.0;
  570.   if (ciCore.tim < 0.0) {
  571.     ciCore.tim += 24.0;
  572.     ciCore.day--;
  573.   }
  574.   CastChart(fTrue);
  575.   mc2 = RFromD(planet[oMC]); k = RFromD(planetalt[oMC]);
  576.   EclToEqu(&mc2, &k);
  577.   cp2 = cp0;
  578.   for (i = 1; i <= cObj; i++) {
  579.     rgalt2[i] = planetalt[i];
  580.   }
  581.  
  582.   /* Loop through 36 hours, dividing it into a certain number of segments. */
  583.   /* For each segment we get the planet positions at its endpoints.        */
  584.  
  585.   for (div = 1; div <= division; div++) {
  586.     ciCore = ciMain;
  587.     ciCore.tim = DecToDeg(ciCore.tim) - 18.0 + 36.0*(real)div/(real)division;
  588.     if (ciCore.tim < 0.0) {
  589.       ciCore.tim += 24.0;
  590.       ciCore.day--;
  591.     } else if (ciCore.tim >= 24.0) {
  592.       ciCore.tim -= 24.0;
  593.       ciCore.day++;
  594.     }
  595.     ciCore.tim = DegToDec(ciCore.tim);
  596.     CastChart(fTrue);
  597.     mc1 = mc2;
  598.     mc2 = RFromD(planet[oMC]); k = RFromD(planetalt[oMC]);
  599.     EclToEqu(&mc2, &k);
  600.     cp1 = cp2; cp2 = cp0;
  601.     for (i = 1; i <= cObj; i++) {
  602.       rgalt1[i] = rgalt2[i]; rgalt2[i] = planetalt[i];
  603.     }
  604.  
  605.     /* During our segment, check to see if each planet rises or sets. */
  606.  
  607.     for (i = 0; i <= cObj; i++) if (!FIgnore(i) && FThing(i)) {
  608.       EclToHorizon(&azi1, &alt1, cp1.obj[i], rgalt1[i], lon, lat, mc1);
  609.       EclToHorizon(&azi2, &alt2, cp2.obj[i], rgalt2[i], lon, lat, mc2);
  610.       j = 0;
  611.       if ((alt1 > 0.0) != (alt2 > 0.0)) {
  612.         d = RAbs(alt1)/(RAbs(alt1)+RAbs(alt2));
  613.         k = Mod(azi1 + d*MinDifference(azi1, azi2));
  614.         j = 1 + (MinDistance(k, rDegHalf) < rDegQuad);
  615.       }
  616.       if (j && occurcount < MAXINDAY) {
  617.         source[occurcount] = i;
  618.         type[occurcount] = j;
  619.         time[occurcount] = 36.0*((real)(div-1)+d)/(real)division*60.0;
  620.         occurcount++;
  621.       }
  622.     }
  623.   }
  624.  
  625.   /* Sort each event in order of time when it happens during the day. */
  626.  
  627.   for (i = 1; i < occurcount; i++) {
  628.     j = i-1;
  629.     while (j >= 0 && time[j] > time[j+1]) {
  630.       SwapN(source[j], source[j+1]);
  631.       SwapN(type[j], type[j+1]);
  632.       SwapR(&time[j], &time[j+1]);
  633.       j--;
  634.     }
  635.   }
  636.  
  637.   /* Now fill out the planet array with the appropriate sector location. */
  638.  
  639.   for (i = 1; i <= cObj; i++) if (!ignore[i] && FThing(i)) {
  640.     planet[i] = 0.0;
  641.     /* Search for the first rising or setting event of our planet. */
  642.     for (s2 = 0; s2 < occurcount && source[s2] != i; s2++)
  643.       ;
  644.     if (s2 == occurcount)
  645.       {
  646. LFail:
  647.       /* If we failed to find a rising/setting bracket around our time, */
  648.       /* automatically restrict that planet so it doesn't show up.      */
  649.       ignore[i] = fTrue;
  650.       continue;
  651.       }
  652. LRetry:
  653.     /* One rising or setting event was found. Now search for the next one. */
  654.     s1 = s2;
  655.     for (s2 = s1 + 1; s2 < occurcount && source[s2] != i; s2++)
  656.       ;
  657.     if (s2 == occurcount)
  658.       goto LFail;
  659.     /* Reject the two events if either (1) they're both the same, i.e. both */
  660.     /* rising or both setting, or (2) they don't bracket the chart's time.  */
  661.     if (type[s2] == type[s1] || time[s1] > 18.0*60.0 || time[s2] < 18.0*60.0)
  662.       goto LRetry;
  663.     /* Cool, we've found our rising/setting bracket. The sector position is */
  664.     /* the proportion the chart time is between the two event times.        */
  665.     planet[i] = (18.0*60.0 - time[s1])/(time[s2] - time[s1])*rDegHalf;
  666.     if (type[s1] == 2)
  667.       planet[i] += rDegHalf;
  668.     planet[i] = Mod(rDegMax - planet[i]);
  669.   }
  670.  
  671.   /* Restore original chart info as we've overwritten it. */
  672.  
  673.   ciCore = ciMain;
  674.   us.fSidereal = fT;
  675. }
  676.  
  677.  
  678. /*
  679. ******************************************************************************
  680. ** Aspect Calculations.
  681. ******************************************************************************
  682. */
  683.  
  684. /* Set up the aspect/midpoint grid. Allocate memory for this array, if not */
  685. /* already done. Allocation is only done once, first time this is called.  */
  686.  
  687. bool FEnsureGrid()
  688. {
  689.   if (grid != NULL)
  690.     return fTrue;
  691.   grid = (GridInfo FPTR *)PAllocate(sizeof(GridInfo), fFalse, "grid");
  692.   return grid != NULL;
  693. }
  694.  
  695.  
  696. /* Indicate whether some aspect between two objects should be shown. */
  697.  
  698. bool FAcceptAspect(obj1, asp, obj2)
  699. int obj1, asp, obj2;
  700. {
  701.   int fSupp;
  702.  
  703.   if (ignorea(asp))    /* If the aspect restricted, reject immediately. */
  704.     return fFalse;
  705.   if (us.fSmartCusp) {
  706.  
  707.     /* Allow only conjunctions to the minor house cusps. */
  708.  
  709.     if ((FMinor(obj1) || FMinor(obj2)) && asp > aCon)
  710.       return fFalse;
  711.  
  712.     /* Prevent any simultaneous aspects to opposing angle cusps,     */
  713.     /* e.g. if conjunct one, don't be opposite the other; if trine   */
  714.     /* one, don't sextile the other; don't square both at once, etc. */
  715.  
  716.     fSupp = (asp == aOpp || asp == aSex || asp == aSSx || asp == aSes);
  717.     if ((FAngle(obj1) || FAngle(obj2)) &&
  718.       (fSupp || (asp == aSqu &&
  719.       (obj1 == oDes || obj2 == oDes || obj1 == oNad || obj2 == oNad))))
  720.       return fFalse;
  721.  
  722.     /* Prevent any simultaneous aspects to the North and South Node. */
  723.  
  724.     if (fSouthNode) {
  725.       if (((obj1 == oNod || obj2 == oNod) && fSupp) ||
  726.         ((obj1 == oSou || obj2 == oSou) && (fSupp || asp == aSqu)))
  727.         return fFalse;
  728.     }
  729.   }
  730.   return fTrue;
  731. }
  732.  
  733.  
  734. /* This is a subprocedure of FCreateGrid() and FCreateGridRelation().   */
  735. /* Given two planets, determine what aspect, if any, is present between */
  736. /* them, and save the aspect name and orb in the specified grid cell.   */
  737.  
  738. void GetAspect(planet1, planet2, ret1, ret2, i, j)
  739. real *planet1, *planet2, *ret1, *ret2;
  740. int i, j;
  741. {
  742.   int k;
  743.   real l, m;
  744.  
  745.   grid->v[i][j] = grid->n[i][j] = 0;
  746.   l = MinDistance(planet2[i], planet1[j]);
  747.   for (k = us.nAsp; k >= 1; k--) {
  748.     if (!FAcceptAspect(i, k, j))
  749.       continue;
  750.     m = l-rAspAngle[k];
  751.     if (RAbs(m) < GetOrb(i, j, k)) {
  752.       grid->n[i][j] = k;
  753.  
  754.       /* If -ga switch in effect, then change the sign of the orb to    */
  755.       /* correspond to whether the aspect is applying or separating.    */
  756.       /* To do this, we check the velocity vectors to see if the        */
  757.       /* planets are moving toward, away, or are overtaking each other. */
  758.  
  759.       if (us.fAppSep)
  760.         m = RSgn2(ret1[j]-ret2[i])*
  761.           RSgn2(MinDifference(planet2[i], planet1[j]))*RSgn2(m)*RAbs(m);
  762.       grid->v[i][j] = (int)(m*60.0);
  763.     }
  764.   }
  765. }
  766.  
  767.  
  768. /* Very similar to GetAspect(), this determines if there is a parallel or */
  769. /* contraparallel aspect between the given two planets, and stores the    */
  770. /* result as above. The settings and orbs for conjunction are used for    */
  771. /* parallel and those for opposition are used for contraparallel.         */
  772.  
  773. void GetParallel(planet1, planet2, planetalt1, planetalt2, i, j)
  774. real *planet1, *planet2, *planetalt1, *planetalt2;
  775. int i, j;
  776. {
  777.   int k;
  778.   real l, alt1, alt2;
  779.  
  780.   l = RFromD(planet1[j]); alt1 = RFromD(planetalt1[j]);
  781.   EclToEqu(&l, &alt1); alt1 = DFromR(alt1);
  782.   l = RFromD(planet2[i]); alt2 = RFromD(planetalt2[i]);
  783.   EclToEqu(&l, &alt2); alt2 = DFromR(alt2);
  784.   grid->v[i][j] = grid->n[i][j] = 0;
  785.   for (k = Min(us.nAsp, aOpp); k >= 1; k--) {
  786.     if (!FAcceptAspect(i, k, j))
  787.       continue;
  788.     l = RAbs(k == aCon ? alt1 - alt2 : RAbs(alt1) - RAbs(alt2));
  789.     if (l < GetOrb(i, j, k)) {
  790.       grid->n[i][j] = k;
  791.       grid->v[i][j] = (int)(l*60.0);
  792.     }
  793.   }
  794. }
  795.  
  796.  
  797. /* Fill in the aspect grid based on the aspects taking place among the */
  798. /* planets in the present chart. Also fill in the midpoint grid.       */
  799.  
  800. bool FCreateGrid(fFlip)
  801. bool fFlip;
  802. {
  803.   int i, j, k;
  804.   real l;
  805.  
  806.   if (!FEnsureGrid())
  807.     return fFalse;
  808.   for (j = 0; j <= cObj; j++) if (!FIgnore(j))
  809.     for (i = 0; i <= cObj; i++) if (!FIgnore(i))
  810.  
  811.       /* The parameter 'flip' determines what half of the grid is filled in */
  812.       /* with the aspects and what half is filled in with the midpoints.    */
  813.  
  814.       if (fFlip ? i > j : i < j) {
  815.         if (us.fParallel)
  816.           GetParallel(planet, planet, planetalt, planetalt, i, j);
  817.         else
  818.           GetAspect(planet, planet, ret, ret, i, j);
  819.       } else if (fFlip ? i < j : i > j) {
  820.         l = Mod(Midpoint(planet[i], planet[j])); k = (int)l;  /* Calculate */
  821.         grid->n[i][j] = k/30+1;                               /* midpoint. */
  822.         grid->v[i][j] = (int)((l-(real)(k/30)*30.0)*60.0);
  823.       } else {
  824.         grid->n[i][j] = SFromZ(planet[j]);
  825.         grid->v[i][j] = (int)(planet[j]-(real)(grid->n[i][j]-1)*30.0);
  826.       }
  827.   return fTrue;
  828. }
  829.  
  830.  
  831. /* This is similar to the previous function; however, this time fill in the */
  832. /* grid based on the aspects (or midpoints if 'acc' set) taking place among */
  833. /* the planets in two different charts, as in the -g -r0 combination.       */
  834.  
  835. bool FCreateGridRelation(fMidpoint)
  836. bool fMidpoint;
  837. {
  838.   int i, j, k;
  839.   real l;
  840.  
  841.   if (!FEnsureGrid())
  842.     return fFalse;
  843.   for (j = 0; j <= cObj; j++) if (!FIgnore(j) || !FIgnore2(j))
  844.     for (i = 0; i <= cObj; i++) if (!FIgnore(i) || !FIgnore2(i))
  845.       if (!fMidpoint) {
  846.         if (us.fParallel)
  847.           GetParallel(cp1.obj, cp2.obj, cp1.alt, cp2.alt, i, j);
  848.         else
  849.           GetAspect(cp1.obj, cp2.obj, cp1.dir, cp2.dir, i, j);
  850.       } else {
  851.         l = Mod(Midpoint(cp2.obj[i], cp1.obj[j])); k = (int)l; /* Calculate */
  852.         grid->n[i][j] = k/30+1;                                /* midpoint. */
  853.         grid->v[i][j] = (int)((l-(real)(k/30)*30.0)*60.0);
  854.       }
  855.   return fTrue;
  856. }
  857.  
  858.  
  859. /*
  860. ******************************************************************************
  861. ** Other Calculations.
  862. ******************************************************************************
  863. */
  864.  
  865. /* Fill out tables based on the number of unrestricted planets in signs by  */
  866. /* element, signs by mode, as well as other values such as the number of    */
  867. /* objects in yang vs. yin signs, in various house hemispheres (north/south */
  868. /* and east/west), and the number in first six signs vs. second six signs.  */
  869. /* This is used by the -v chart listing and the sidebar in graphics charts. */
  870.  
  871. void CreateElemTable(pet)
  872. ET *pet;
  873. {
  874.   int i, s;
  875.  
  876.   ClearB((lpbyte)pet, (int)sizeof(ET));
  877.   for (i = 0; i <= cObj; i++) if (!FIgnore(i)) {
  878.     pet->coSum++;
  879.     s = SFromZ(planet[i]);
  880.     pet->coSign[s-1]++;
  881.     pet->coElemMode[(s-1)&3][(s-1)%3]++;
  882.     pet->coElem[(s-1)&3]++; pet->coMode[(s-1)%3]++;
  883.     pet->coYang += (s & 1);
  884.     pet->coLearn += (s < sLib);
  885.     if (!FCusp(i)) {
  886.       pet->coHemi++;
  887.       s = inhouse[i];
  888.       pet->coHouse[s-1]++;
  889.       pet->coModeH[(s-1)%3]++;
  890.       pet->coMC += (s >= sLib);
  891.       pet->coAsc += (s < sCan || s >= sCap);
  892.     }
  893.   }
  894.   pet->coYin   = pet->coSum  - pet->coYang;
  895.   pet->coShare = pet->coSum  - pet->coLearn;
  896.   pet->coDes   = pet->coHemi - pet->coAsc;
  897.   pet->coIC    = pet->coHemi - pet->coMC;
  898. }
  899.  
  900. /* calc.c */
  901.